This is an R Markdown Notebook. When you execute code within the notebook, the results appear beneath the code.
Try executing this chunk by clicking the Run button within the chunk or by placing your cursor inside it and pressing Ctrl+Shift+Enter.
library(qiime2R)
Registered S3 method overwritten by 'data.table':
method from
print.data.table
Registered S3 method overwritten by 'htmlwidgets':
method from
print.htmlwidget tools:rstudio
library(tidyr)
library(microbiome)
Loading required package: phyloseq
Loading required package: ggplot2
microbiome R package (microbiome.github.com)
Copyright (C) 2011-2019 Leo Lahti,
Sudarshan Shetty et al. <microbiome.github.io>
Attaching package: ‘microbiome’
The following object is masked from ‘package:ggplot2’:
alpha
The following object is masked from ‘package:base’:
transform
library(compositions)
Welcome to compositions, a package for compositional data analysis.
Find an intro with "? compositions"
Attaching package: ‘compositions’
The following objects are masked from ‘package:stats’:
anova, cor, cov, dist, var
The following objects are masked from ‘package:base’:
%*%, norm, scale, scale.default
library(Dict)
library(reshape2)
Attaching package: ‘reshape2’
The following object is masked from ‘package:tidyr’:
smiths
library(heatmaply)
Loading required package: plotly
Attaching package: ‘plotly’
The following object is masked from ‘package:ggplot2’:
last_plot
The following object is masked from ‘package:stats’:
filter
The following object is masked from ‘package:graphics’:
layout
Loading required package: viridis
Loading required package: viridisLite
Registered S3 method overwritten by 'dendextend':
method from
rev.hclust vegan
Registered S3 method overwritten by 'seriation':
method from
reorder.hclust vegan
======================
Welcome to heatmaply version 1.3.0
Type citation('heatmaply') for how to cite the package.
Type ?heatmaply for the main documentation.
The github page is: https://github.com/talgalili/heatmaply/
Please submit your suggestions and bug-reports at: https://github.com/talgalili/heatmaply/issues
You may ask questions at stackoverflow, use the r and heatmaply tags:
https://stackoverflow.com/questions/tagged/heatmaply
======================
Attaching package: ‘heatmaply’
The following object is masked from ‘package:compositions’:
normalize
library(KEGGREST)
#This code block is to stage all the files needed for visualization
#Stage taxonomy file using Qiime2R
taxa_vsearch <- parse_taxonomy(taxonomy = read_qza('../taxonomy_vsearch/classification.qza')$data)
taxa_blast <- parse_taxonomy(taxonomy = read_qza('../taxonomy_blast/classification-blast.qza')$data)
#Stage metadata file
sample_metadata <- read.csv('../read-files/metadata.tsv', sep='\t')
#Stage unstratified/strat kegg pathway files
unstrat_kegg <- read.csv('../picrust2-strat/out/KEGG_pathways_out/path_abun_unstrat.tsv', sep = '\t')
strat_kegg <- read.csv('../picrust2-strat/out/KEGG_pathways_out/path_abun_strat.tsv', sep = '\t')
#Stage unstrat/strat KEGG ortholog files
ko_unstrat <- read.csv('../picrust2-strat/out/KO_metagenome_out/unstrat_ko_metagenome.tsv', sep = '\t')
ko_strat <- read.csv('../picrust2-strat/out/KO_metagenome_out/strat_ko_metagenome.tsv', sep = '\t')
#Open selected KOs (for plotting relevant Kegg pathways)
deg_kos <- read.csv('./select_kos.txt', sep = '\t', header = TRUE)
rownames(deg_kos) <- deg_kos[,1]
deg_kos[,1] <- NULL
#open the pathway to ko file
pth_2_ko <- read.csv('~/bioinfo_pipelines/picrust2-2.4.2/picrust2/default_files/pathway_mapfiles/KEGG_pathways_to_KO.tsv', header =FALSE, sep='\t')
#Filter to only Select KOs
pth_2_ko <- pth_2_ko[pth_2_ko$V1 %in% rownames(deg_kos),]
#Retrieve Sample ids so I can call it as a variable
sample_ids <- c(as.character(sample_metadata$Sample.id))
#retrieve BRITE hierarchy from KEGG using KEGGREST
keggGet(dbentries = 'br08901')
Error in .getUrl(url, .flatFileParser) : Not Found (HTTP 404).
#This code block is to load the required files to plot NSTI vs relative abundance
abs_count_otu <- read_qza(file = '../dada2/table.qza')$data
rel_count_otu <- as.data.frame(apply(abs_count_otu,2,function(x){x/sum(x)}))
nsti_asvs <- as.data.frame(read.csv('../picrust2-strat/out/marker_predicted_and_nsti.tsv.gz', sep='\t', row.names = 1))
rel_nsti <- merge(nsti_asvs,rel_count_otu,by="row.names",all.x=TRUE)
rel_nsti <- melt(rel_nsti, measure.vars = sample_ids, variable.name = 'sample.id', value.name = 'rel.abundance')
#Group accding to consortia
ggplot_rel_plot <- ggplot(data=rel_nsti,mapping= aes(x=metadata_NSTI, y=rel.abundance, color=sample.id)) + geom_point()
plot(ggplot_rel_plot)